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Abstract 



The short-time critical dynamics of propagation of 
damage in the Ising ferromagnet in two dimensions is 
studied by means of Monte Carlo simulations. Start- 
ing with equilibrium configurations at T = oo and 
magnetization M — 0, an initial damage is created 
by flipping a small amount of spins in one of the 
two replicas studied. In this way, the initial dam- 
age is proportional to the initial magnetization Mq 
in one of the configurations upon quenching the sys- 
tem at Tc, the Onsager critical temperature of the 
ferromagnetic-paramagnetic transition. It is found 
that, at short times, the damage increases with an 
exponent 9d = 1.915(3), which is much larger than 
the exponent 9 — 0.197 characteristic of the initial in- 
crease of the magnetization M{t). Also, an epidemic 
study was performed. It is found that the average dis- 
tance from the origin of the epidemic {{E?{t))) grows 
with an exponent z* Ki -q p^ 1.9, which is the same, 
within error bars, as the exponent 9d- However, 
the survival probability of the epidemics reaches a 
plateau so that (5 = 0. On the other hand, by quench- 
ing the system to lower temperatures one observes 
the critical spreading of the damage at To — O.SlTp, 
where all the measured observables exhibit power 
laws with exponents 9d — 1.026(3), 5 = 0.133(1), 
and z* = 1.74(3). 
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1 Introduction 

The critical behavior of a statistical system near a 
second-order phase transition is characterized by crit- 
ical exponents defining a universality class [1]. This 
behavior is due to the fact that the characteristic 
length of the physical system given by the spatial 
correlation length (^) becomes infinite at Tc, the crit- 
ical point, as (T — Tc)~'^ , where v is the correlation 
length exponent, leading to the observation of scale 
invariance. But at T = Tc, the correlation length 
grows according to ^(t) oc i^^^, where z is the dy- 
namic critical exponent, so for a finite system of side 
L and large enough time, ^ saturates at a certain cor- 
relation time r, such that t (x L^. Because z > 0, 
the correlation time increases rapidly with the lattice 
size, and therefore, in Monte Carlo simulations it is 
very difficult to generate independent configurations. 
This effect, known as critical slowing down, affects 
not only measurements in equilibrium but also mea- 
surements of dynamic observables. In this context, 
the development of the short-time dynamic (STD) 
theory [2] provides a very useful tool not only to avoid 
the critical slowing down but also for the measure- 
ments of critical exponents. In fact, 20 years ago, 
using renormalization group techniques, Jansenn et 
al. [2] showed that the scaling behavior is not only 
valid in equilibrium but also that dynamic scaling 
holds during the short-time regime of the evolution 
of a critical system. It is worth mentioning that re- 
cent progress in the understanding of phase transi- 
tions and critical phenomena has been boosted by a 
large number of investigations based on the applica- 



tion of STD. Indeed, the STD approach has proved 
to be an invaluable tool for the study of a large vari- 
ety of equilibrium systems (e.g., models such as the 
Ising Ong, the XY [UMi], the Heisenberg [H], the 
Potts [TT1I271I25] . the Blume-Capel pilSU] . etc.). Fur- 
thermore, the generalization of the STD scaling ap- 
proach to far-from-equilibrium |31| and self-organized 
critical [32J systems poses a large theoretical chal- 
lenge. 

In order to apply the STD scaling approach in 
Monte Carlo simulations of a magnetic system, the 
sample should be prepared at a very high tempera- 
ture, with a small magnetization Mq remaining. It 
is then suddenly quenched to the critical tempera- 
ture Tc and released to the selected dynamic evolu- 
tion of the model (Glauber, Metropolis, heat bath, 
etc.). The dynamic scaling relationship obtained by 
Jansenn et al. [2!, in the case of systems of finite 
size, and for the fcth moment of the magnetization, 
reads |MI 



Summing up, after a microscopic time such as 
C(^mic) is of the order of the lattice spacing, one ob- 
serves an initial increase of the magnetization accord- 
ing to Eq. (21) up to a mesoscopic time of the order 
while for t > tm the standard relax- 
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ation of the order parameter, namely, M 
observed. 

On the other hand, the damage spreading (DS) 
method is a standard technique frequently used to 
study the propagation of perturbations in lattice sys- 
tems. In past decades it has been applied to very dif- 
ferent systems such as the Ising magnet, spin glasses, 
Potts models, cellular automata, biological models, 
etc. (for more details of this technique, see the recent 
review |53| and references therein). The methodol- 
ogy used to study the propagation of a perturbation 
by means of DS is to start the simulations with two 
well-equilibrated configurations, S^i^^ and S-^{T), 
which differ from one another only in the state of a 
small number of sites. Then, both configurations are 
allowed to evolve in time, with the same sequence 
of random numbers (in terms of the simulations this 



mC^) (t, e, L, Mo) = b ^^'"m''^'^ (6 H, b^'^e, b ^L, 5^« Afc^idition means that both configurations experience 



(1) 
where e = 1 — T/Tq is the reduced critical temper- 
ature, and /3 and v are the equilibrium critical ex- 
ponents for the order parameter and the correlation 
length, respectively; z is the dynamic critical expo- 
nent and the new independent exponent xq is the 
scaling dimension of the initial magnetization Mq. 

Numerical simulations support the theoretical pre- 
diction for the short-time dynamic scaling |33fl40) . 
Further, the short-time dynamic scaling is found to 
be very general [MIH1[171[31I1TH51] . 

In the thermodynamic limit (L — > oo) and choosing 
the scaling factor b = t^^^ (so that the first argument 
of the scaling function on the right-hand side of Eq. 
^ is set to 1), the magnetization of the system {k = 
1) can be written as 



the same thermal noise). Thus, at certain time t, the 
difference between S^ and S^ will only be as a con- 
sequence of the small initial perturbation or damage 
introduced in the system at time t = 0. 

In order to measure the difference between these 
two configurations, the Hamming distance or damage 
(D(i)) is defined [54J as 
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where the exponent 6 = (xq — /3/i')/z has been intro- 
duced. In most cases 6* > 0, i.e., for e = the small 
initial magnetization increases in the short-time re- 
gion. 



where the summation runs over the total number of 
sites N, and the index I {1 < I < N) is the label that 
identifies the sites of the configurations. S'^itjT) is 
an equilibrium configuration of the system at temper- 
ature T and time t, while S^{t,T) is the perturbed 
configuration [54,'55J. 

If the initial perturbation introduced in S^ is small 
{D{t = 0) — > 0), at least two possible scenarios are ex- 
pected: (i) D{t — ^ oo) goes to a finite nonzero value, 
and the perturbation is relevant to the system, or (ii) 
the small perturbation vanishes after a certain time 
and Z)(i — )■ oo) — > 0. This behavior introduces a new 



continuous and irreversible critical transition between 
a state where damage heals and a state where the 
perturbation propagates into the system. This tran- 
sition is known as damage transition and, in the case 
of the Ising model in two dimensions with Glauber 
dynamics, it occurs at a damage critical temperature 
Td ^ 0.992(2)Tc [53], where Tc is the Onsager crit- 
ical temperature for the paramagnetic-ferromagnetic 
transition. 

The universality class of the damage transition is 
still an open question. Grassberger [ST] conjectures 
that the DS transition may belong to the directed 
percolation (DP) universality class [58H60] if its crit- 
ical point does not coincide with a critical transi- 
tion of the physical system, e.g., the critical temper- 
ature of the Ising magnet. Numerical simulations of 
different systems such as the Domany-Kinzel cellu- 
lar automata, the two-dimensional (2D) Ising model 
with Swendsen-Wang dynamics [51], as well as the 
deterministic cellular automata with small noise |62| 
provide evidence that the DS transition is charac- 
terized by critical exponents of the DP universality 
class. However, there are also other systems where 
the DS transition has a non-DP behavior, such as 
the case of the KaufFman model [63, 64' or the Ziff- 
Gulari-Barshad (ZGB) model in 2D [65J (for more 
details, see |53)). 

STD has been applied to DS in the Ising model at 
the Onsager temperature with heat bath dynamics 
in two and three dimensions by Grassberger [39] . He 
found that the damage exhibits power-law behavior 
of the form 

0(1)0^"°, (4) 

where 6o is the initial increase exponent for the dam- 
age. The exponent Od was found to be 9d = 0.191(3) 
{d = 2) and 9d = 0.104(3) (d = 3). Also, the dy- 
namic dependence of the survival probability of dam- 
age P{t), defined as the probability at time t of find- 
ing nonzero damage, was analyzed in [391. Grass- 
berger found, as expected [66], a power-law decrease 
with an exponent S, of the form 
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the spins damaged on average) a decrease in damage 
is observed with an exponent 6'jj = 0.43(2) (d = 3). 

Within this context, the aim of this work is to re- 
port results obtained upon the study of the STD of 
damage spreading in the two-dimensional Ising model 
with Metropolis dynamics. In contrast to previous 
studies |53) where the initial configuration is equili- 
brated at a given temperature where the DS is sub- 
sequently measured, here we equilibrate the system 
at T = oo (i.e., we obtain fully uncorrelated initial 
configurations), and afterwards we quench the system 
at the temperature of interest for the measurement 
of the damage. Therefore, we expect to gain insight 
into the propagation of the damage in a quite dif- 
ferent scenario than in standard studies, and search 
for a relationship between damage spreading and the 
short-time dynamic behavior of relevant observables, 
e.g. the magnetization. 

The paper is organized as follows. In Sec. [51 we 
briefly describe the Ising Model and give the simu- 
lation details, in Sec. [3] we report and discuss the 
results obtained, and finally in Sec. JU the conclu- 
sions are presented. 

2 Ising Model 

As was mentioned in the previous section, we study 
the STD of damage spreading in the two-dimensional 
classical Ising model ^7\. For this reason, in the 
present section we introduce a brief description of 
this archetypical system largely used to study phase 
transitions and critical phenomena in magnetic sys- 
tems. 

In the Ising model, each site ai of the lattice of 
size L X L represents a spin variable, which interacts 
with its nearest-neighbor spins with a constant of ex- 
change J. If J > 0, and in the absence of an external 
magnetic field, the Hamiltonian ("H) can be written 
as 
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(6) 



with 6 « 0.9 {d = 2) and S w 1.1 {d = 3). For ran- 
dom independent initial configurations (i.e., half of 



where ai is the Ising spin variable that can assume 
two different values ai — ±1, the indexes I < i,j < N 
are used to label the spins, and the summation runs 



over all the nearest-neighbor pairs of spins. In the 
absence of an external magnetic field and at low tem- 
perature, the system is, for d > 1, in the ferromag- 
netic phase and, on average, the majority of spins 
are pointing in the same direction. In contrast, at 
high temperature the system maximizes the entropy, 
thermal fluctuations break the order and the system 
is in the paramagnetic phase. This ferromagnetic- 
paramagnetic critical transition is a second-order 
phase transition and it occurs at a well-defined crit- 
ical temperature (Tc). In the two-dimensional case, 
one has exactly kTc/J = 2/ln{l + V2) = 2.269..., 
where k is the Boltzmann constant. 

With the aim of studying the STD of DS, we start 
the simulation with a configuration S with strictly 
zero initial magnetization. This constraint was in- 
troduced by starting the simulation with all spins 
of configuration S"^ equal to cr = f , and then N/2 
spins (randomly chosen) of the lattice are flipped to 
cr = — 1. In this context, the initial configuration 
with uncorrelated spins and M = corresponds to 
T — oo. Att = 0, a replica S^ was created, and the 
initial damage was introduced by flipping Nq spins 
randomly chosen. In this way, the initial damage of 
the configuration S^ is related to the nonzero initial 
magnetization as Mq — 2D{t = 0) = 2Nq/L'^. 

In order to set the time scale, we assume that dur- 
ing a Monte Carlo time step (mcs) all spins of the 
system (L x L in total) have the chance to be flipped 
once, on average. 

In this context, the spreading of damage is stud- 
ied by assuming Metropolis dynamics and the simu- 
lations were performed on the square lattice of size 
L X L (100 < L < 4096), with periodic boundary 
conditions between the borders of the lattice. 



3 Results and discussion 

In order to test the computational protocol, in Fig. 
[T]we show the results obtained for the initial increase 
in the magnetization as a function of time (Eq. ©), 
for lattices of different size L, T = Tc, and using 
Metropolis algorithm. In this example, Mq = 0.01 is 
assumed. The best fit of the data yields 9 = 0.196(1) 
in excellent agreement with the best available value 
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Figure 1: Time evolution of magnetization at the 
critical temperature Tc for the Ising model by ap- 
plying Metropolis algorithm. Results were obtained 
for Mq = 10~^ and different lattice sides L. The 
dashed line corresponds to the best fit of the data 
and has slope 6 = 0.196(1). 



for the Ising model-with Metropolis dynamic-in the 
Mo -^ limit, namely, 6 = 0.197(1) [33J. 

On the other hand. Fig. [UJa) shows plots of 
D{t) versus t obtained by keeping T — Tc con- 
stant and varying the initial magnetization Mq (and 
therefore the initial damage D{0)), for a lattice of 
size L = 1024. The dashed line corresponds to a 
power-law behavior (Eq. S]), but the results obtained 
show that the value of the exponent 9d depends on 
the initial magnetization (Mq). The obtained val- 
ues of 0£) also depend on the lattice size (not shown 
here for the sake of space), e.g., for a lattice side 
L = 1024 and T = Tc, 6*13(^0 = lO'^) = 0.81(1), 
6»£,(Mo = 10-3) = 1.21(1) and OoiMo = 9 x 10"^) = 
1.84(1) (see Figured^a)), and for L ^ 512, eniMo = 
10-2) = 0.76(1), 0d{Mq = 10-3) = 1.15(1) and 
6»£,(Afo = 9 X 10-6) = 1.77(1). For this reason and 
in order to establish the value of this exponent in 
the limits Mq — >■ and L — ^ 00, we develop an epi- 
demic study by flipping only the five central spins of 
the lattice. So, the initial magnetization for the lat- 
tice of size L X L is Mq = 2D{t = 0) = 2 x 5/L'^. 



The results obtained are shown in Fig. [5] (b) for 
different lattice sizes as indicated in figure and at 
the critical temperature Tq- The inset of Fig. [Hb) 
shows the scaled results. After a transient period of 
t ~ 10^ mcs, the damage exhibits a power-law behav- 
ior, but the value of the exponent 9d also depends 
on the lattice size L, as is shown in Fig. [3] The 
data obtained suggest a size dependence of the form 
9d{L) = 9d{L — ^ oo) 4- aL^^, and the best fit of the 
data gives OoiL -j> cxd) = 1.915(3) (see Fig. E]). As 
can be seen this exponent is much bigger (roughly 
one order of magnitude) than the exponent for the 
initial increase of the magnetization, — 0.197 |33| . 

In order to establish the dependence of this expo- 
nent on the dynamic rules applied to the lattice, we 
also performed simulations of the short-time dynam- 
ics of damage with heat bath, Glauber and Metropo- 
lis dynamics, in the epidemic case. The results ob- 
tained are shown in Fig. 21 As is expected, the dam- 
age heals for heat bath dynamics at T = Tq, and it 
spreads in the cases of Metropolis and Glauber dy- 
namics. In both cases, the value of the exponent is 
given by 9d = 1.91(1), in complete agreement with 
our previous result. 

On the other hand, we analyze the dynamic depen- 
dence of the damage survival probability {P{t)). Fig- 
ure [5] shows the results obtained for different temper- 
atures and L = 2048. As can be observed, P{t) — > 
for T -^ Tc, while P(t) reaches a stationary value 
at Tc ■ This behavior suggests that the damage tran- 
sition occurs at a temperature lower than Tc, and 
at Tc the decay exponent (see Eq. ([S])) is given by 
S — 0. This result is in agreement with that reported 
by Montani et al. [68] in the two-dimensional Ising 
model with Glauber dynamics. 

On the other hand, we also studied the average dis- 
tance from the origin of the propagation of damage, 
defined as 139 



(R^t)) = 



EAisHt,T)-s^{t,T)y\ 

j:i{iS^it,T)-S^{t,T))) 



(7) 



where r = \{x — L/2)\. 

Figure [5] shows a log-log plot of (i?^(t)) versus t, 
for lattices of different size L and T = Tc- These 
data suggest a power-law behavior for (i?^(i)) of the 
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Figure 2: (Color online) (a) Log-log plot of damage vs 
t, obtained for a lattice of side L ~ 1024, by applying 
the Metropolis algorithm and keeping T = Tc con- 
stant. Different values of the initial magnetization 
AIq = 2D{0) are used, as listed in the main panel. 
The dashed line has slope 9d{L) = 1.84(1). (b) Dam- 
age vs t for the epidemic study, with Mq — 2 x 5/i^, 
T ~ Tc and different lattice sizes, as shown in the 
figure. Inset: Collapse of the data. 
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Figure 3: Linear behavior of the exponent O^, as a 
function of L~^ , the inverse of the lattice size. The 
asymptotic hmit values 9d{L — ^ oo) — 1.915(3). 
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Figure 5: (Color online) Log-log plot of the survival 
probability of the damage {P{t)) vs t obtained for a 
lattice size L — 2048, and different temperatures T, 
as indicated in the figure. As can be observed, at T = 
Tc, P{t) reaches a stationary value corresponding to 
5 = 0. The fuU line {T = To ^ O.SlTc) has slope 
(5 = 0.133(1). 
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Figure 4: (Color online) Log-log plot of damage as 
a function of time, for a lattice of size L = 4096, 
T = Tc, Mo = ^(O) = 5/L^ and different dynamic 
rules applied, as indicated in the figure. The dashed 
line has slope 6'£) sa 1.91. 
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where the best fit of the data (in the limit L —>■ cx)) 
gives z* = 1.92(4). This exponent is similar, within 
error bars, to the exponent Od. 

In order to understand the spatiotemporal evo- 
lution of both the magnetization and the damage, 
we also analyzed the profiles along the x direction. 
For this purpose, we studied a strip geometry of size 
L X M with strictly zero magnetization. The initial 
damage was introduced by flipping only three spins 
of the central column (i = Af/2) in the configuration 
labeled S'5(t, Tc), so Mq = 2D{0) = 2 x 3/(L x M). 
Thus, we define the damage and magnetization pro- 
files as 
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Figure 6: (Color online) Average distance from the 
origin of the propagation of damage R^{t) vs t, for 
Tc and different values of the lattice size L. The 
dashed line has slope z = 1.92. 

where Sfj{t,T) and Sfj{t,T) are the reference and 
damaged configurations at site I of coordinates {i, j}, 
respectively. By assuming these definitions, the pro- 
files M{i,t) and D{i,t) represents the average mag- 
netization and damage of the zth (i = 1, ...,M) row 
of the system, which runs parallel to the x direction, 
respectively. 

Figure [Tlja) shows the results obtained for a lattice 
of size L X M = 100 X 500, at T = T^ and different 
times as indicated in the figure caption. As can be 
observed, the propagation of damage is faster than 
that of the magnetization. This result is in complete 
agreement with the fact that for the initial increase 
exponents one has 9d > 0. Also, in Fig. [lib) we 
show, in more detail, the profiles for the initial times 
of the propagation, namely, t < 100 mcs. 

As was mentioned in Sec. [TJ after a mesoscopic 
time tm ~ Mq^ (y = z/xq), the magnetization de- 
cays as M ^ t~^/'^^ (see Fig. [5]). For the Ising model 
y K. 3.96 and P/vz « 0.0625, so that for an initial 
magnetization A/q = 0.1, i„j ^ lO'^mcs, Mq = 0.05, 
tm ~ lO^mcs, and Mq — 10^^, t^ ^ lO^^mcs. On 
the other hand, after a certain time t*, the damage 
reaches a saturation value. This characteristic time 
t* was determined as a function of lattice size and the 
initial magnetization Afo (see the inset of Fig. El) . As 
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Figure 7: Damage and magnetization profiles at the 
X direction for a lattice of size L x M = 100 x 
500. The results are obtained using Metropolis al- 
gorithm, at T = Tc and for different times: (a) 
t = 1, 100, 400, 1024, 2025, and 2500 mcs; (b) the ini- 
tial times t = 1,4,9,16,25, and 36 mcs. The full 
line indicates, in both cases, the initial damage (or 
magnetization) . 
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Figure 8: (Color online) Time evolution of the mag- 
netization at the critical temperature Tc for the Ising 
model with Mq = 0.05 and 0.1. The dashed line 
corresponds to the initial increase M ~ i^, with 
6 ss 0.19, and the full line shows the standard re- 
laxation of M(t) after a mesoscopic time t„j, namely 
M - t-/^/"'', with 13/ vz = 0.0625. The inset 
shows the crossover time of damage (i*) and mag- 
netization (tm) vs A/q. The dotted line has slope 
—z/xq w —3.96, while the slope of the dash-dotted 
line is -0.53(1). 



can be seen, there are no finite size effects (at least for 
the lattice sizes used in the simulations). Moreover, 
the inset of Fig. [8] shows a power-law behavior of the 
form t* ^ Mq^", with yn = 0.53(1). By comparing 
this value with the magnetization given by y w 3.96 
(squares in the inset of Fig. [5]), it can be observed 
that the damage spreads over the lattice much more 
quickly than the magnetization, in agreement with 
previous results shown in this paper. 

Finally, the behavior observed for the survival 
probability (see Fig. [5]) suggests that the critical tem- 
perature for damage transition is lower than Tc- In 
fact, there is another (lower) temperature at which 
both D{t) and P{t) behave as power laws. This hap- 
pens at Td = 0.51(l)Tc and the corresponding ex- 
ponents are Od — 1.026(3) and S — 0.133(1) (see 
Figs. [5]and[9]). Also, the exponent z* = 1.74(3) was 
determined at To- 



Figure 9: (Color online) Log- log plot of the damage 
D{t) vs t obtained for a lattice size L = 2048 and 
different temperatures T as indicated in the figure. 
The dashed line (corresponding to T = Tc) has slope 
OoiL = 2048) = 1.91, and the fuh line {T = Td ^ 
0.51Tc) has slope 0d = 1.026(3). 

4 Conclusions 

The propagation of damage and the short-time crit- 
ical dynamics of the d = 2 Ising model are studied. 
For this purpose, fully uncorrelated initial configura- 
tions (T = oo) are suddenly quenched at the mea- 
sured temperature. The proposed methodology for 
the study of the spreading of the damage is in con- 
trast to previous studies [53] where the equilibration 
and measurement temperatures are identical. It is 
found that by quenching the system at Tc, or close 
to it, the spreading of the damage is much faster 
than the initial increase of the magnetization. How- 
ever, the critical temperature for damage spreading 
is far below Tc, namely close to To = 0.51 (l)Tc. 
It is worth mentioning that the values of the ex- 
ponents measured at To, namely, 9u = 1.026(3), 

5 = 0.133(1), and z* = 1.74(3), are different from 
those reported by Grassberger for the Ising model 
with heat bath dynamics at Tc [39]. However, it must 

be mentioned that both systems are quite different. 
In our case, the damage was introduced in a com- 
pletely uncorrelated configuration that was suddenly 
quenched at a temperature T. Subsequently, the 



damage starts to evolve with the standard Metropo- 
Us dynamics at T. Therefore, a completely different 
scenario is expected: the critical behavior is not ob- 
served at Tc-or close to it, as usual-but at a much 
lower temperature given by T ~ O.SlTc- This result 
would mean that the damage spreading Transition 
observed in our case is completely different, and for 
this reason, the values of the critical exponents re- 
ported here are not the same as those reported by 
Grassberger |39) . 

Finally, we would like to remark that our results 
may indicate that the damage spreading critical tem- 
perature would systematically decrease when increas- 
ing the equilibration temperature. To = 0.51 (l)Tc 
being the lower bound as obtained for TequU = oo. 
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